---
title: "CARRS RDD analysis using a standardized cutoff"
knit: (function(input, ...) {rmarkdown::render(input, output_file=".pdf")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r , setup, include=FALSE, eval = TRUE}

knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""

)
 knitr::opts_knit$set(
   root.dir = ""    
 )
```

```{r packages}
rm(list=ls())
  library(tidyverse)
  library(gridExtra)
  library(dplyr)
  library(rdrobust)
  library(rddensity)
  library(viridisLite)
  library(foreign)
  library(ggplot2)
  library(data.table)
  library(doBy)
  library(socviz)
  library(openxlsx) 
  library(modelsummary)
  library(labelled)
  library(gtsummary)
  library(gt)
  library(kableExtra)
  library(ggpubr)
```

```{r loaddata}
data_treatment<-data.frame(read.csv("carrs_full.csv"))
data_diagnosis<-data.frame(read.csv("carrs_selected.csv"))
```

```{r fuc}
##Defining the common function##
  carrs_rdd = function(a,b,d,outcome,od,yl,boxh,boxv,clr) {
    
    d$Y<-a
    d$X<-b
    gname=od
    yl=yl*100
    runningvar=b
    out<- rdrobust(d$Y, d$X, c = 0, 
                    cluster =d$pid, 
                    covs = cbind(d$w01, d$w23, d$w45), 
                    kernel = "triangular", 
                    p = 1,rho=1)

    ret <- data.frame(
    sex = c(formatC(round(100*out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(100*out$ci[3,1],2),2,format="f"), " - ", formatC(round(100*out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2])))
    
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
    output<<-ret %>% add_rownames 

    d$bw <- out$bws[1,1]
    bw = out$bws[1,1]
    
   plotsub<-d[runningvar>=-bw & runningvar<=bw,] ##the plot is only for the values within the bandwidth
   
   c=0
  
  plot1 =rdplot(plotsub$Y, plotsub$X,
                  kernel = "triangular", p=1, 
                  covs = cbind(plotsub$w01, plotsub$w23, plotsub$w45), 
                         ci=95, hide=TRUE)
  rdplot_mean_bin = plot1$vars_bins[,"rdplot_mean_bin"]
  rdplot_mean_y   = plot1$vars_bins[,"rdplot_mean_y"]
  y_hat           = plot1$vars_poly[,"rdplot_y"]
  x_plot          = plot1$vars_poly[,"rdplot_x"]

  rdplot_mean_y = rdplot_mean_y*100
  y_hat = y_hat*100

  y_hat_r=y_hat[x_plot>=c & x_plot!=0]
  y_hat_l=y_hat[x_plot<c & x_plot!=0]
  x_plot_r=x_plot[x_plot>=c & x_plot!=0]
  x_plot_l=x_plot[x_plot<c & x_plot<0]

  x.label="\nStandardized blood pressure\n"
  y.label=outcome
  x.lim=c(min(plotsub$X, na.rm=T),max(plotsub$X, na.rm=T))
  y.lim=c(0,15)
  
  color_plot<- ggplot() + 
    geom_point(aes(x = rdplot_mean_bin, y = rdplot_mean_y, 
                    ), size=2,  colour ="#999999", na.rm = TRUE) +
    geom_line(aes(x = x_plot_l, y = y_hat_l,colour ='x_plot_l'),
                col= clr, size=1,na.rm = TRUE) +
    geom_line(aes(x = x_plot_r, y = y_hat_r,colour ='x_plot_r'), 
                col= clr, size=1, na.rm = TRUE) +

    labs(x = x.label, y = y.label) +
    labs(caption = paste("MSE optimal bandwidth: +/- ", round(bw,2),"\n")) +
    labs(y = y.label, x = x.label) +
    theme(legend.position = "None") +
    geom_vline(xintercept = c, size = 0.5) +
    coord_cartesian(xlim = x.lim, ylim = y.lim) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"),
          axis.text=element_text(family="serif", size=28), 
          axis.title=element_text(family="serif", size=28),
          axis.title.x = element_text(vjust = 5),
          axis.title.y = element_text(vjust = -5),
          plot.caption = element_text(family="serif", size=24, vjust = 10))

    out<<-out
    bw<<-bw 
    color_plot<<- color_plot
  }
```

**Graphs

\newpage

**Outcome : Diagnosis**
**Females and males**  

```{r dfull}
data_diagnosis<-as.data.table(data_diagnosis)

carrs_rdd(data_diagnosis$hypdiag,data_diagnosis$runningvar,data_diagnosis,"Probability (%)\n\n","overall_diagnosis.pdf",0.1,-0.9,10.5,"#25AC82FF")
diag_total <-output%>%mutate(Total=sex)%>%select(-sex)
bwTdiag  <- bw 
plotdiagT<-color_plot
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```
\newpage

**Outcome : Diagnosis**

**Female**  

```{r dfemale}
data_diagnosis1<-data_diagnosis[female==1,]

carrs_rdd(data_diagnosis1$hypdiag,data_diagnosis1$runningvar_female,data_diagnosis1,"Probability (%)\n\n", "female_diagnosis.pdf",0.1,-0.85,10.5,"#46337EFF")
diag_female <-output%>%mutate(Female=sex)%>%select(-sex)
bwFdiag  <- bw
plotdiagF<-color_plot
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```
\newpage

**Outcome : Diagnosis (Males)**

```{r dmale}
data_diagnosis2<-data_diagnosis[female==0,]

carrs_rdd(data_diagnosis2$hypdiag,data_diagnosis2$runningvar_male,data_diagnosis2,"Probability (%)\n\n","male_diagnosis.pdf",0.1,-0.93,10.5,"#B63679FF")
diag_male <-output%>%mutate(Male=sex)%>%select(-sex)
bwMdiag  <- bw
plotdiagM<-color_plot
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```

```{r ,results= TRUE, fig.height=22}
figurediag <- ggarrange(plotdiagT, plotdiagF, plotdiagM,
                    labels = c("Total population", "Females", "Males"),
                    ncol = 1,
                    font.label = list(size=28, face="bold", family="serif"),
                    label.x = c(-0.12,-0.03,-0.01), label.y = c(1.1,1.1,1.1)) +
  theme(plot.margin = margin(2,1,2,1, "cm")) 

figurediag
```



\newpage

**Outcome : Treatment**

**Females and males**  

```{r tfull}
data_treatment<-as.data.table(data_treatment)

carrs_rdd(data_treatment$hypdrug,data_treatment$runningvar,data_treatment,"Probability (%)\n\n","overall_treatment.pdf",0.16,-0.83,16.8,"#25AC82FF")
treat_total <-output%>%mutate(Total=sex)%>%select(-sex)
plottreatT<-color_plot
bwTtreat  <- bw
rm(bw)
```

```{r ,results= TRUE}
summary(out)
```

\newpage

**Outcome : Treatment**

**Female**  

```{R tfemale}
data_treatment1<-data_treatment[female==1,]

carrs_rdd(data_treatment1$hypdrug,data_treatment1$runningvar_female,data_treatment1,"Probability (%)\n\n", "female_treatment.pdf",0.16,-0.91,16.8,"#46337EFF")
treat_female <-output%>%mutate(Female=sex)%>%select(-sex)
plottreatF<-color_plot
bwFtreat  <- bw
rm(bw)
```



```{r ,results= TRUE}
summary(out)
```

\newpage

**Outcome : Treatment (Male)**

```{r tmale}
data_treatment2<-data_treatment[female==0,]

carrs_rdd(data_treatment2$hypdrug,data_treatment2$runningvar_male,data_treatment2,"Probability (%)\n\n","male_treatment.pdf",0.16,-0.95,16.8,"#B63679FF")
treat_male <-output%>%mutate(Male=sex)%>%select(-sex)
plottreatM<-color_plot
bwMtreat  <- bw 
rm(bw)
```


```{r ,results= TRUE}
summary(out)
```

```{r ,results= TRUE, fig.height=22}
figuretreat <- ggarrange(plottreatT, plottreatF, plottreatM,
                    labels = c("Total population", "Females", "Males"),
                    ncol = 1,
                    font.label = list(size=28, face="bold", family="serif"),
                    label.x = c(-0.12,-0.03,-0.01), label.y = c(1.1,1.1,1.1)) +
                    theme(plot.margin = margin(2,1,2,1, "cm")) 

figuretreat
```


```{r ,results= TRUE, fig.height=25,fig.width=20}
figurecomb <- ggarrange(figurediag, figuretreat,
                    labels = c("Diagnosis", "Treatment"),
                    ncol = 2,
                    font.label = list(size=32, face="bold", family="serif"),
                    label.x = c(0.4,0.37), label.y = c(1.03,1.03)) +
                    theme(plot.margin = margin(2,1,2,1, "cm")) 
figurecomb
filename <- paste0("combined.pdf")
    ggsave(filename, dpi = 300)
    
filename <- paste0("combined.jpeg")
    ggsave(filename, dpi = 300)
```


```{r}
total_table<-merge(diag_total,treat_total, by="rowname")
total_table<-total_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))


colnames(total_table)[2] <- "Diagnosis"   
colnames(total_table)[3] <- "Treatment"   

female_table<-merge(diag_female,treat_female, by="rowname")
female_table<-female_table%>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(female_table)[2] <- "Diagnosis"   
colnames(female_table)[3] <- "Treatment"    

male_table<-merge(diag_male,treat_male, by="rowname")
male_table<-male_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(male_table)[2] <- "Diagnosis" 
colnames(male_table)[3] <- "Treatment"    
```
##table making function


```{r}
#combining diagnosis and treatment
  filename <- paste0("mainresults_combined_updated.tex")
  caption=paste0("Effect of home-based screening on hypertension diagnosis and treatment")

  comb_table<-rbind(total_table,female_table,male_table)

combined<-kbl(comb_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccc", format = "latex",col.names = c(" ","Diagnosis","Treatment"),digits = 2) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

combined

```

##Tables

\newpage

```{r ,results= FALSE}
tables = function(tdata) {

set_gtsummary_theme(list(
  "tbl_summary-fn:percent_fun" = scales::label_percent(accuracy = 0.1)
))
  
tout<<-  tdata %>%  mutate(total = 1) %>% select(total, everything())  %>% select(-female) %>% 
  select(-c('runningvar_female','runningvar_male','runningvar')) %>%
  tbl_summary(
  statistic = list(all_continuous() ~ "{mean} ({sd})",
                   all_categorical() ~ "{p}%",
                    total ~ "{n}"),
  digits = list(all_categorical() ~ c(1,1),
                all_continuous() ~ c(0,0)),
  missing = "no" ,
  label = list(total="Number",
            age_calc="Mean age, years (SD)",
            agecatf="Age groups",
            religion="Religion",
            eduf="Education categories",
            diabetesknow="Suffers from diabetes",
            heartknow="Suffered a heart attack ",
            strokeknow="Suffered a stroke"))  

}
```

##Diagnosis

```{r ,results= FALSE}
#Hypertension diagnosis
tables.data<-data_diagnosis[,c("age_calc","agecatf","religion","eduf","diabetesknow","heartknow", "strokeknow","female","runningvar","runningvar_female","runningvar_male")] 
```

```{r ,results= FALSE}
#overall diagnosis 
tdata<- tables.data
tables(tdata)
overall_diag<-tout %>%
  modify_header(update = stat_0 ~ "**Overall**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```

```{r ,results= FALSE}
#female diagnosis 
tdata<- tables.data %>% filter(female==1) 
tables(tdata)
fem_overall_diag<-tout %>%
  modify_header(update = stat_0 ~ "**Overall**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```

```{r ,results= FALSE}
#male diagnosis 
tdata<-tables.data %>% filter(female==0)
tables(tdata)
male_overall_diag<-tout %>%
  modify_header(update = stat_0 ~ "**Overall**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```


```{r ,results= FALSE}
#analytical diagnosis 
tdata<- tables.data  %>%  filter(runningvar>=-bwTdiag & runningvar<=bwTdiag)  
tables(tdata)
analytical_diag<-tout%>%
  modify_header(update = stat_0 ~ "**Within BW**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)

```

```{r ,results= FALSE}
#analytical female diagnosis 
tdata<-tables.data  %>% filter(female==1)  %>% filter(runningvar>=-bwFdiag & runningvar<=bwFdiag)  
tables(tdata)
fem_analytical_diag<-tout%>%
  modify_header(update = stat_0 ~ "**Within BW**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```

```{r ,results= FALSE}
#analytical male diagnosis 
tdata<- tables.data  %>% filter(female==0)  %>% filter(runningvar>=-bwMdiag & runningvar<=bwMdiag)  
tables(tdata)
male_analytical_diag<-tout %>%
  modify_header(update = stat_0 ~ "**Within BW**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```

##Treatment tables

```{r ,results= FALSE}
#Hypertension treatment

tables.data<-data_treatment[,c("age_calc","agecatf","religion","eduf","diabetesknow","heartknow", "strokeknow","female","runningvar","runningvar_female","runningvar_male")] 
```

```{r ,results= FALSE}
#overall treatment 
tdata<- tables.data 
tables(tdata)
overall_treat<-tout %>%
  modify_header(update = stat_0 ~ "**Overall**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```

```{r ,results= FALSE}
#female treatment 


tdata<- tables.data %>% filter(female==1) 
tables(tdata)
fem_overall_treat<-tout %>%
  modify_header(update = stat_0 ~ "**Overall**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)


```

```{r ,results= FALSE}
#male treatment 
tdata<-tables.data %>% filter(female==0) 
tables(tdata)
male_overall_treat<-tout %>%
  modify_header(update = stat_0 ~ "**Overall**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```

```{r ,results= FALSE}
#analytical treatment 

tdata<- tables.data %>% filter(runningvar>=-bwTtreat & runningvar<=bwTtreat)  
tables(tdata)
analytical_treat<-tout%>%
  modify_header(update = stat_0 ~ "**Within BW**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA)
rm(tout,tdata)
```

```{r ,results= FALSE}
#analytical female treatment 

tdata<-tables.data  %>% filter(female==1)  %>% filter(runningvar>=-bwFtreat & runningvar<=bwFtreat)  
tables(tdata)
fem_analytical_treat<-tout%>%
  modify_header(update = stat_0 ~ "**Within BW**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA)
rm(tout,tdata)
```

```{r ,results = FALSE}
#analytical male treatment 

tdata<- tables.data  %>% filter(female==0) %>%  filter(runningvar>=-bwMtreat & runningvar<=bwMtreat)  
tables(tdata)
male_analytical_treat<-tout %>%
  modify_header(update = stat_0 ~ "**Within BW**") %>% 
  modify_footnote(c(all_stat_cols()) ~ NA) 
rm(tout,tdata)
```


```{r ,results= FALSE}
# table with male, female and total population
table_final_diag_all<-tbl_merge(
  tbls = list(overall_diag, analytical_diag, fem_overall_diag,fem_analytical_diag,male_overall_diag,male_analytical_diag))%>%
  modify_spanning_header(
    c(stat_0_1 , stat_0_2) ~ "**Total**", 
    c(stat_0_3 , stat_0_4) ~ "**Female**", 
    c(stat_0_5 , stat_0_6) ~ "**Male**")

as_kable_extra(table_final_diag_all, format = "latex",booktabs = T,  escape = FALSE,
  linesep = '', caption = "Characteristics of the overall and within optimal bandwidth samples, outcome: hypertension diagnosis")  %>%
  kable_styling(position = "center",latex_options = "HOLD_position")%>%
   footnote(general=c("The table displays the characteristics of the sample used in the analysis with hypertension diagnosis as the outcome. The overall sample consists of all individuals meeting the eligibility criteria. The Within BW sample consists of those individuals with a standardized blood pressure within the mean squared error optimal bandwidth.", "Abbreviation: BW = bandwidth, SD = standard deviation."), general_title="",threeparttable = T, fixed_small_size = T) %>%
  save_kable(file = "sumchar_diagnosis_all.tex",format="latex")

```


```{r ,results = FALSE}
table_final_treat_all<-tbl_merge(
  tbls = list(overall_treat,analytical_treat,fem_overall_treat,fem_analytical_treat,male_overall_treat,male_analytical_treat)) %>%
  modify_spanning_header(
    c(stat_0_1 , stat_0_2) ~ "**Total**", 
    c(stat_0_3 , stat_0_4) ~ "**Female**", 
    c(stat_0_5 , stat_0_6) ~ "**Male**")

as_kable_extra(table_final_treat_all, format = "latex",booktabs = T, escape = FALSE,
  linesep = '', caption = "Characteristics of the overall and within optimal bandwidth samples, outcome: hypertension treatment")   %>%
  kable_styling(position = "center",latex_options = "HOLD_position")%>%
   footnote(general=c("The table displays the characteristics of the sample used in the analysis with hypertension treatment as the outcome. The overall sample consists of all individuals meeting the eligibility criteria. The Within BW sample consists of those individuals with a standardized blood pressure within the mean squared error optimal bandwidth.", "Abbreviation: BW = bandwidth, SD = standard deviation."),general_title="",threeparttable = T, fixed_small_size = T) %>%
  save_kable(file = "sumchar_treatment_all.tex",format="latex",escape = FALSE)

```
**Summary characteristics - Hypertension diagnosis**
```{r ,results = 'asis'}
table_final_diag_all%>%
  as_gt()
```

**Summary characteristics - Hypertension treatment**

```{r ,results = 'asis'}
table_final_treat_all%>%
  as_gt()
```

``` {r RDsensitivityBWfunction2, eval = TRUE}
 bwsens <- function(Yi,bw) {
    out =  rdrobust(y = Yi, X, kernel = "triangular", p = 1, h = bw,  all = T,
                    c = 0, cluster = clus, covs = cov, rho = 1 )
    summary(out)
    
    ret <- data.frame(
      res = c(as.character(round(100*out$coef[1,1],2)), as.character(round(out$pv[3,1],2)), 
              paste0( "(", as.character(round(100*out$ci[3,1],2)), " - ", as.character(round(100*out$ci[3,2],2)), ")" ), 
              as.character(round(out$bws[1,1],2)), as.character(out$N_h[1] + out$N_h[2]))
    )
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
    ret
  
    return(list(ret=ret))
  }
```

``` {r SensistivityBwDiag, eval = TRUE }
## Total population
  X = data_diagnosis$runningvar
  Y = data_diagnosis$hypdiag
  clus = data_diagnosis$pid
  cov = cbind(data_diagnosis$w01, data_diagnosis$w23, data_diagnosis$w45)
  
  # lists with bandwidth +/-10% of optimal bandwidth in steps of 0.05
  bbwwT = rep(seq(from=bwTdiag-0.2, to=bwTdiag+0.2, by=0.05), 1)
  bbwwT

  Tdiag <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(Tdiag) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(bbwwT)) {
  assign("Tdiagres", bwsens(Yi = Y,bw=bbwwT[i]))
  Tdiag[i] <- Tdiagres
  Tdiag[i] <- as.matrix(Tdiag[i])
  }

## Men
  X = data_diagnosis2$runningvar_male
  Y = data_diagnosis2$hypdiag
  clus = data_diagnosis2$pid
  cov = cbind(data_diagnosis2$w01, data_diagnosis2$w23, data_diagnosis2$w45)
  
  # lists with bandwidth +/-10% of optimal bandwidth in steps of 0.05
  bbwwM = rep(seq(from=bwMdiag-0.2, to=bwMdiag+0.2, by=0.05), 1)
  bbwwM

  Mdiag <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(Mdiag) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(bbwwM)) {
  assign("Mdiagres", bwsens(Yi = Y,bw=bbwwM[i]))
  Mdiag[i] <- Mdiagres
  Mdiag[i] <- as.matrix(Mdiag[i])
  }
  
## Women
  X = data_diagnosis1$runningvar_female
  Y = data_diagnosis1$hypdiag
  clus = data_diagnosis1$pid
  cov = cbind(data_diagnosis1$w01, data_diagnosis1$w23, data_diagnosis1$w45)
  
  # lists with bandwidth +/-10% of optimal bandwidth in steps of 0.05
  bbwwF = rep(seq(from=bwFdiag-0.2, to=bwFdiag+0.2, by=0.05), 1)
  bbwwF

  Fdiag <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(Fdiag) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(bbwwF)) {
  assign("Fdiagres", bwsens(Yi = Y,bw=bbwwF[i]))
  Fdiag[i] <- Fdiagres
  Fdiag[i] <- as.matrix(Fdiag[i])
  }
```

``` {r SensistivityBwTreat, eval = TRUE }
## Full sample
  X = data_treatment$runningvar
  Y = data_treatment$hypdrug
  clus = data_treatment$pid
  cov = cbind(data_treatment$w01, data_treatment$w23, data_treatment$w45)
  
  # lists with bandwidth +/-10% of optimal bandwidth in steps of 0.05
  bbwwT = rep(seq(from=bwTtreat-0.2, to=bwTtreat+0.2, by=0.05), 1)
  bbwwT

  Ttreat <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(Ttreat) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(bbwwT)) {
  assign("Ttreatres", bwsens(Yi = Y,bw=bbwwT[i]))
  Ttreat[i] <- Ttreatres
  }
  
## Men
  X = data_treatment2$runningvar_male
  Y = data_treatment2$hypdrug
  clus = data_treatment2$pid
  cov = cbind(data_treatment2$w01, data_treatment2$w23, data_treatment2$w45)
  
  # lists with bandwidth +/-10% of optimal bandwidth in steps of 0.05
  bbwwM = rep(seq(from=bwMtreat-0.2, to=bwMtreat+0.2, by=0.05), 1)
  bbwwM

  Mtreat <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(Mtreat) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(bbwwM)) {
  assign("Mtreatres", bwsens(Yi = Y,bw=bbwwM[i]))
  Mtreat[i] <- Mtreatres
  }
  
  ## Women
  X = data_treatment1$runningvar_female
  Y = data_treatment1$hypdrug
  clus = data_treatment1$pid
  cov = cbind(data_treatment1$w01, data_treatment1$w23, data_treatment1$w45)
  
  # lists with bandwidth +/-10% of optimal bandwidth in steps of 0.05
  bbwwF = rep(seq(from=bwFtreat-0.2, to=bwFtreat+0.2, by=0.05), 1)
  bbwwF

  Ftreat <- data.frame(matrix(, nrow=5, ncol=0))
  rownames(Ftreat) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  for (i in 1:length(bbwwF)) {
  assign("Ftreatres", bwsens(Yi = Y,bw=bbwwF[i]))
  Ftreat[i] <- Ftreatres
  }
``` 

```{r }
##COMBINED - DIAGNOSIS
  filename <- paste0("bw_diag_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension diagnosis with varying bandwidths")

comb_table<-rbind(as.matrix(Tdiag),as.matrix(Fdiag),as.matrix(Mdiag))

combined<-kbl(comb_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccccccc", format = "latex",digits = 2, col.names = c("(1)", "(2)", "(3)", "(4)","(5)","(6)","(7)","(8)","(9)")) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10)  %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15)  %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level. The bandwidth was varied by +/- 0.2 in intervals of 0.05, column (5) contains the results."," Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

```



```{r }
##COMBINED- TREATMENT
  filename <- paste0("bw_treat_combined.tex")
  caption=paste0("Effect of home-based screening on hypertension treatment with varying bandwidths")

comb_table<-rbind(as.matrix(Ttreat),as.matrix(Ftreat),as.matrix(Mtreat))

combined<-kbl(comb_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccccccc", format = "latex",digits = 2, col.names = c("(1)", "(2)", "(3)", "(4)","(5)","(6)","(7)","(8)","(9)")) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Total", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Female", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10)  %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Male", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15)  %>%
    footnote(general = c("The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and CIs resulting from robust bias-corrected standard errors clustered at the individual-level. The bandwidth was varied by +/- 0.2 in intervals of 0.05, column (5) contains the results.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

combined

```


``` {r PlottingBPvalues, eval = TRUE}
sysplot <- ggplot(data = data_treatment, aes(sys)) +
    geom_histogram(binwidth = 1, fill = "#7f7f7f") +

    geom_vline(xintercept=139.5, color="#228D8DFF", lwd=1.2) + 
    scale_x_continuous(breaks=seq(70,230,10), expand = c(0.01,0)) + # ticks of x-axis
    labs(x = "\n\n", y = "Frequency") + 
    scale_y_continuous(breaks=seq(0,600,100), expand = c(0,0)) + 

  # remove grid and grey background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text=element_text(family="serif", size=34), axis.title=element_text(family="serif", size=36)) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

  sysplot

diasplot <- ggplot(data = data_treatment, aes(dias)) +
    geom_histogram(binwidth = 1, fill = "#7f7f7f") +

    geom_vline(xintercept=89.5, color="#228D8DFF", lwd=1.2) + 
    scale_x_continuous(breaks=seq(40,150,10), expand = c(0.01,0)) + # ticks of x-axis
    labs(x = "", y = "Frequency") +
    scale_y_continuous(breaks=seq(0,800,100), expand = c(0,0)) + 

  # remove grid and grey background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text=element_text(family="serif", size=34), axis.title=element_text(family="serif", size=36)) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

  diasplot
```

```{r ,results= TRUE, fig.height=25,fig.width=20}
sysdias <- ggarrange(sysplot, diasplot,
                    labels = c("Systolic blood pressure", "Diastolic blood pressure"),
                    nrow = 2,
                    font.label = list(size=36, face="bold", family="serif"),
                    label.x = c(-0.07,-0.08), label.y = c(1.10,1.10)) +
                    theme(plot.margin = margin(3,1,2,1, "cm")) 

sysdias
filename <- paste0("sysdias.pdf")
    ggsave(filename, dpi = 300)
filename <- paste0("sysdias.jpeg")
    ggsave(filename, dpi = 300)
```

``` {r densitytests, eval = TRUE}
trtdens <-rddensity(data_treatment$runningvar, c=0, p=2, q=3, all=TRUE)
summary(trtdens)

trtdensplot <- rdplotdensity(trtdens,data_treatment$runningvar,histFillCol = 4,histBreaks=seq(-3,3,0.1))
trtdensplot 

filename <- paste0("DensityTrt.jpeg")
    ggsave(filename, width = 5, height = 4, dpi = 300, units = "in")

    
diagdens <-rddensity(data_diagnosis$runningvar, c=0, p=2, q=3)
summary(diagdens)

diagdensplot <- rdplotdensity(diagdens,data_diagnosis$runningvar,histFillCol = 4,histBreaks=seq(-5,5,0.1))
diagdensplot 

filename <- paste0("DensityDiag.jpeg")
    ggsave(filename, width = 5, height = 4, dpi = 300, units = "in")

```



